Homicides in Santa Catarina: a Random Forest Application
by William Suzuki and Ian Fukushima
The objective of this article is to explore and model the dataset of homicides in the state of Santa Catarina. After request, the state’s security bureau of Santa Catarina gave us the data. Ian Fukushima, my colleague, inspired me to make the maps and analysis of this data, I’m grateful for his support.
Load, Explore and Transform
Let’s start by loading packages:
#preliminaries -----------------------------------------
library(knitr) #for kable, stylized table in r markdown
library(tidyverse) #for pipe
library(dplyr) #select() function
library(sp)
library(spdep)
library(rgdal)
library(Amelia) #missmap()
library(ggmap)
library(rgeos) #gUnaryUnion
library(xts) #extended time series
library(magick) # to make gifs
library(randomForest)
library(gstat) #idw function
library(tmap) #to use tmap function
library(randomForest)
library(grid) #to plot two tmaps together
library(raster) #raster()Include some paths to save images. You must change those paths for your machine
path1 = "/home/william/Dropbox/working/projects_git/Crime-SC/tests william/"
path_root = "/home/william/Dropbox/working/projects_git/Crime-SC/"
path_data = "/home/william/Dropbox/working/projects_git/Crime-SC/Dados"load.csv() will load the data. Ian used Google Maps’ API to find the geographic coordinates. The raw data that came only with the address of each homicide incident.
#prepare base----------------------------------------
#load base
crime = read.csv(file=paste0(path_data,"/df_sc_crime_geocoded.csv"), header=TRUE, sep=",", fileEncoding="latin1")crime %>% class## [1] "data.frame"The data have 10745 incidents of homicide.
crime %>% dim()## [1] 10745 25The following are the variables. We are going to use only "Data", "lat","long".
crime %>% names()## [1] "X" "Crime"
## [3] "Data" "Hora"
## [5] "ID" "Endereço"
## [7] "Município" "Bairro"
## [9] "Tipo.Logradouro" "Logradouro"
## [11] "Número" "Tipo.de.local"
## [13] "Localidade" "PosiTime"
## [15] "lat" "long"
## [17] "accuracy" "formatted_address"
## [19] "address_type" "status"
## [21] "index" "MOTIVAÇÃO"
## [23] "VIOLÊNCIA.DOMÉSTICA" "MENOSP.DISCRIM.CONDIÇÃO.DE.MULHER"
## [25] "ESTUPRO..ANTECEDENTE"Here is the summary of our data. We are going to use "Data", "lat","long". But variables as "Municipio","Número","Endereço","Bairro" are important when using Google’s API for geocoding.
summary(crime)## X Crime
## Min. : 1 HOMICÍDIO :9162
## 1st Qu.: 2687 LATROCÍNIO : 665
## Median : 5373 PESSOA MORTA POR POLICIAL MILITAR EM SERVIÇO : 619
## Mean : 5373 LESÃO CORPORAL SEGUIDA DE MORTE : 170
## 3rd Qu.: 8059 PESSOA MORTA POR POLICIAL CIVIL EM SERVIÇO : 63
## Max. :10745 PESSOA MORTA POR POLICIAL MILITAR FORA DE SERVIÇO: 50
## (Other) : 16
## Data Hora ID
## 6/12/2015: 12 21:00 : 638 Min. : 1
## 13/3/2010: 11 23:00 : 604 1st Qu.: 2687
## 13/5/2018: 11 22:00 : 579 Median : 5373
## 31/3/2018: 11 20:00 : 530 Mean : 5373
## 1/1/2013 : 10 0:00 : 512 3rd Qu.: 8059
## 11/7/2010: 10 19:00 : 438 Max. :10745
## (Other) :10680 (Other):7444
## Endereço
## , JOINVILLE, Santa Catarina, BRASIL : 50
## , FLORIANÓPOLIS, Santa Catarina, BRASIL : 22
## RUA SEIS DE JANEIRO , JOINVILLE, Santa Catarina, BRASIL: 21
## , SÃO JOSÉ, Santa Catarina, BRASIL : 15
## , ITAJAÍ, Santa Catarina, BRASIL : 13
## , CRICIÚMA, Santa Catarina, BRASIL : 12
## (Other) :10612
## Município Bairro Tipo.Logradouro
## JOINVILLE :1133 CENTRO : 990 RUA :7234
## FLORIANÓPOLIS:1063 INTERIOR : 321 ESTRADA :1105
## ITAJAÍ : 561 : 216 RODOVIA : 701
## CHAPECÓ : 487 RURAL : 193 AVENIDA : 642
## SÃO JOSÉ : 467 MONTE ALEGRE : 161 : 311
## CRICIÚMA : 380 PARANAGUAMIRIM: 142 SERVIDÃO: 282
## (Other) :6654 (Other) :8722 (Other) : 470
## Logradouro Número Tipo.de.local
## : 322 :7686 VIA PÚBLICA :5728
## GERAL : 185 20 : 21 RESIDÊNCIA :2956
## BR 101 : 71 200 : 20 OUTRO : 848
## BR 470 : 35 30 : 15 ESTABELECIMENTO COMERCIAL: 489
## SANTA CATARINA: 30 55 : 14 BARES E SIMILARES : 452
## GETÚLIO VARGAS: 28 12 : 13 ESTABELECIMENTO PENAL : 122
## (Other) :10074 (Other):2976 (Other) : 150
## Localidade PosiTime
## ZONA URBANA :9851 2014-08-31 03:00:00: 6
## ZONA RURAL : 882 2014-09-06 22:00:00: 6
## NÃO INFORMADA : 4 ZONA URBANA : 6
## \\"CANTINA PER TUTI\\",","VIA PÚBLICA: 1 2014-03-20 19:00:00: 5
## \\"E\\",","VIA PÚBLICA : 1 2015-02-26 04:00:00: 5
## 406, \\"G\\",","VIA PÚBLICA : 1 2016-06-12 06:00:00: 5
## (Other) : 5 (Other) :10712
## lat long accuracy
## -28.6010292: 91 -49.151029 : 91 route :5971
## -26.2834213: 73 -48.8452269: 73 street_address:2136
## -26.6140628: 41 -53.7265491: 41 locality : 790
## -27.6140791: 27 -48.6370861: 27 establishment : 781
## -26.0713285: 25 -48.6177384: 25 premise : 371
## (Other) :10471 (Other) :10471 (Other) : 679
## NA's : 17 NA's : 17 NA's : 17
## formatted_address
## Estrada Geral, Santa Catarina, Brazil : 91
## Joinville - State of Santa Catarina, Brazil : 74
## BR-101, Santa Catarina, Brazil : 41
## Estrada de Linha, Paraíso - SC, 89906-000, Brazil : 41
## São José, Santa Catarina - Barreiros, São José - State of Santa Catarina, Brazil: 27
## (Other) :10454
## NA's : 17
## address_type
## route :5971
## street_address :2136
## locality,political : 790
## premise : 371
## administrative_area_level_2,political: 321
## (Other) :1139
## NA's : 17
## status
## OK :10720
## Error - NA - Verifify : 17
## route : 3
## Chapecó, State of Santa Catarina, Brazil : 1
## clothing_store,establishment,home_goods_store,point_of_interest,store: 1
## establishment,point_of_interest : 1
## (Other) : 2
## index MOTIVAÇÃO VIOLÊNCIA.DOMÉSTICA
## OK : 5 NÃO INFORMADA :3480 NÃO :10178
## 1 : 1 DESAVENÇA :2853 NÃO INFORMADO: 1
## 10 : 1 TRÁFICO DE DROGAS:1726 SIM : 566
## 100 : 1 PASSIONAL : 846
## 1000 : 1 AÇÃO POLICIAL : 745
## 10000 : 1 ROUBO : 655
## (Other):10735 (Other) : 440
## MENOSP.DISCRIM.CONDIÇÃO.DE.MULHER ESTUPRO..ANTECEDENTE
## NÃO :4489 NÃO :10559
## NÃO INFORMADO:6232 NÃO INFORMADO: 160
## SIM : 24 SIM : 26
##
##
##
## kable(crime[1:3,], caption="original data of homicides") #shows stylized table of data frame| X | Crime | Data | Hora | ID | Endereço | Município | Bairro | Tipo.Logradouro | Logradouro | Número | Tipo.de.local | Localidade | PosiTime | lat | long | accuracy | formatted_address | address_type | status | index | MOTIVAÇÃO | VIOLÊNCIA.DOMÉSTICA | MENOSP.DISCRIM.CONDIÇÃO.DE.MULHER | ESTUPRO..ANTECEDENTE |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | HOMICÍDIO | 1/1/2008 | 0:00 | 1 | RUA MONTEVIDÉU , CHAPECÓ, Santa Catarina, BRASIL | CHAPECÓ | FORTES | RUA | MONTEVIDÉU | VIA PÚBLICA | ZONA URBANA | 2008-01-01 00:00:00 | -27.0996274 | -52.5998438 | route | R. Montevidéu, Chapecó - SC, Brazil | route | OK | 1 | TRÁFICO DE DROGAS | NÃO | NÃO INFORMADO | NÃO | |
| 2 | HOMICÍDIO | 1/1/2008 | 0:00 | 2 | RUA PEGASUS 187, JOINVILLE, Santa Catarina, BRASIL | JOINVILLE | JARDIM PARAÍSO | RUA | PEGASUS | 187 | OUTRO | ZONA URBANA | 2008-01-01 00:00:00 | -26.2144603 | -48.8141942 | route | R. Pegasus - Jardim Paraiso, Joinville - SC, 89240-000, Brazil | route | OK | 2 | TRÁFICO DE DROGAS | NÃO | NÃO INFORMADO | NÃO |
| 3 | HOMICÍDIO | 1/1/2008 | 6:00 | 3 | AVENIDA AMÂNDIO CABRAL , BALNEÁRIO BARRA DO SUL, Santa Catarina, BRASIL | BALNEÁRIO BARRA DO SUL | CENTRO | AVENIDA | AMÂNDIO CABRAL | VIA PÚBLICA | ZONA URBANA | 2008-01-01 06:00:00 | -26.4594161 | -48.6030489 | route | Av. Amândio Cabral, Balneário Barra do Sul - SC, 89247-000, Brazil | route | OK | 3 | TRÁFICO DE DROGAS | NÃO | NÃO INFORMADO | NÃO |
Many variables came in factor format. We transform all variables that are factor into character.
#transform all variables that are factor into character
crime2 = crime
len1 = length(names(crime2))
ii = 1
for (ii in 1:len1) {
vv = crime[,ii]
if (is.factor(vv)) {
vv1 <- as.character(vv)
crime2[,ii] = vv1
}
}#X and ID show the same thing
identical(crime2$X,crime2$ID)## [1] TRUE#years and days, and hours of crime
crime2$Data %>% class## [1] "character"#we can aggregate the data into months or trimestersTake a look at the types of homicide. We’ll consider all of them as just homicide.
kable(crime2$Crime %>% table)| . | Freq |
|---|---|
| HOMICÍDIO | 9162 |
| INFANTICÍDIO | 8 |
| LATROCÍNIO | 665 |
| LESÃO CORPORAL SEGUIDA DE MORTE | 170 |
| PESSOA MORTA POR POLICIAL CIVIL EM SERVIÇO | 63 |
| PESSOA MORTA POR POLICIAL CIVIL FORA DE SERVIÇO | 8 |
| PESSOA MORTA POR POLICIAL MILITAR EM SERVIÇO | 619 |
| PESSOA MORTA POR POLICIAL MILITAR FORA DE SERVIÇO | 50 |
There are only some lines that are missing, but for variables that we don’t use.
#map of missing in crime2 from Amelia
missmap(crime2)
#transform lat and long into numeric
crime2$lat=as.numeric(crime2$lat)
crime2$long=as.numeric(crime2$long) #there are a few missings, 25 and 19
is.na(crime2$lat) %>% sum## [1] 25is.na(crime2$long) %>% sum## [1] 19#drop places with missing lat and long
mis_lat = is.na(crime2$lat)
mis_lon = is.na(crime2$lat)
crime3 = crime2[!(mis_lat|mis_lon),]
crime3 %>% dim()## [1] 10720 25#we drop 25 observationsThere are some outliers for latitude, for example, -27027.
summary(crime3$lat)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -27027.00 -27.59 -27.05 -34.48 -26.74 45.14#there are some outliers in lat and long
#proportion of points that does not make sense
summary(crime3$long)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -119.70 -49.58 -48.83 -49.40 -48.65 35.37len2 = dim(crime3)[1]
len2## [1] 10720#10720 number of lines in datasum(crime3$long < -54)## [1] 8#8 number of outliers
sum(crime3$long < -54)/len2## [1] 0.0007462687#0.0007462687
#-54 long is more than the most western frontiers of Santa Catarina
#-53.8376609 most west of SC
#-47.2330741 most east of SCsum(crime3$long > -47)## [1] 63#latitude explore outlierssummary(crime3$lat)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -27027.00 -27.59 -27.05 -34.48 -26.74 45.14#look at google maps to see which values makes sense
# -25.3345494 most north of SC
# -29.5 most south of SCMaps
We use function get_map() from ggmap to get a “screen shot” from Google maps.
# maps ------------------------------------
#google maps
map = get_map(location=c(-53.9, -29.5 ,-47,-25.33),maptype = 'roadmap')
plot(map)
#transform crime3 into spatial points data.frame
coordinates(crime3) = ~long+lat
plot(crime3$long,crime3$lat )
#plot of coordinates doesn't make sense because of outliers
plot(crime3)
Those strange results are important to be shown. Those are some examples of problems that we face using real world data.
crime3 %>% class## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"The projection of a spatial object is important. Below we see that crime3 doesn’t have a projection.
crime3@proj4string #there is no projection## CRS arguments: NAHere we give a projection for crime3
crime3@proj4string = CRS("+proj=longlat +datum=WGS84")
plot(crime3) #there are some outliers in the west, it doesn't make sense
points_occur = data.frame(lon=coordinates(crime3)[,1],
lat=coordinates(crime3)[,2])(crime_sc1=ggmap(map, ylab = "Latitude", xlab = "Longitude") +
#geom_polygon(aes(x=long, y=lat, group=group),
# fill='grey', size=.2, color='blue', data=, alpha=0)+
geom_point(data=points_occur, aes(x=lon, y=lat), size=.5, color='red', alpha=.5)
)
#save image to disk
ggsave("crime_sc1.jpg", width = 20, height = 20, units = "cm",dpi = 600)Municipal Borders
#import data of municipalities:
munic = readOGR("~/Dropbox/working/RAW_DATA/shpbrasil/municip07.shp", layer="municip07")## OGR data source with driver: ESRI Shapefile
## Source: "/home/william/Dropbox/working/RAW_DATA/shpbrasil/municip07.shp", layer: "municip07"
## with 5561 features
## It has 3 fields#pick up only municipalities that are in the state of Santa Catarina
state = substr(munic$CODIGOIB1,1,2)
munic=munic[state=='42',]
munic@proj4string = CRS("+proj=longlat +datum=WGS84")
plot(munic)
#plot map with municipaility borders and points of crime and google maps
munic_fortify = fortify(munic, region="CODIGOIB1")
munic_fortify %>% class## [1] "data.frame"#we first transform spatial polygons into data frame
(crime_sc2=ggmap(map, ylab = "Latitude", xlab = "Longitude") +
geom_polygon(data=munic_fortify, color='black', alpha=0,
aes(y=lat, x=long, group=group, fill='grey'))+
geom_point(data=points_occur, aes(x=lon, y=lat), size=.75, color='red', alpha=.5)
)
#save image
ggsave("crime_sc2.jpg", width = 20, height = 20, units = "cm",dpi = 600)#we need to exclude points outside santa catarina's borders
#get outer border of Santa Catarina
munic_dissolve = gUnaryUnion(munic)
plot(munic_dissolve)
Note that munic_dissolve is a spatial polygon, then we use it to select values of crime3 which is spatial points. So this is a way to crop values in a spatial points.
#exclude outliers, we solve the problem
crime4 = crime3[munic_dissolve,]
plot(crime4)
Now that we excluded all outliers we have a recognizable map.
Hexagon Polygons
#create hexagon polygons --------------------------------------------------
#code from:
#http://strimas.com/spatial/hexagonal-grids/#f1
size = 0.3
set.seed(123)
hex_points = spsample(munic_dissolve, type = "hexagonal", cellsize = size)
hex_points %>% class #spatial points## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"hex_points %>% plot #just make points in a hexagon fashion
hex_points %>% length #there are 109 points ## [1] 109hex_grid = HexPoints2SpatialPolygons(hex_points, dx = size)
hex_grid %>% class #from spatial points to spatial polygons## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"plot(hex_grid)
#number of hexagons in map
len_hex = length(hex_grid) ; len_hex ## [1] 109plot(munic_dissolve, col = "grey50", bg = "light blue", axes = TRUE)
plot(hex_points, col = "black", pch = 20, cex = 0.5, add = T)
plot(hex_grid, border = "orange", add = T)
#to save image
jpeg("hexagon_sc.jpg")
plot(munic_dissolve, col = "grey50", bg = "light blue", axes = TRUE)
plot(hex_points, col = "black", pch = 20, cex = 0.5, add = T)
plot(hex_grid, border = "orange", add = T)
dev.off() ## png
## 2hex_fortify = fortify(hex_grid)
(hex_sc2=ggmap(map, ylab = "Latitude", xlab = "Longitude") +
geom_polygon(data=hex_fortify, color='black', alpha=0,
aes(y=lat, x=long, group=group, fill='grey'))+
geom_point(data=points_occur, aes(x=lon, y=lat), size=.75, color='red', alpha=.5)
)
Note that there are some points that are inside any hexagon. Those points will not be accounted in the final data frame. This is a drawback in using hexagon polygons.
ggsave("hex_sc2.jpg", width = 20, height = 20, units = "cm",dpi = 600)#another view of the data
plot(hex_grid)
plot(crime4, col="red" , add=TRUE)
#prepare spatial polygons data frame
crime4$count = 1
names(crime4)[24] #"count"## [1] "count"sp1 = crime4[,24] #take just the "count" variable
sp1 %>% class## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"res = over(hex_grid, sp1 , fn=sum) #count nº of homicides in each hexagon
res %>% dim## [1] 109 1is.na(res$count) %>% sum## [1] 5res$count[is.na(res$count)] = 0 #hexagons with no crime have 0 crimeshex_sp = SpatialPolygonsDataFrame(hex_grid, data=data.frame(ID=1:length(hex_grid)), match.ID = FALSE)
hex_sp$crime = res$count
#plot of crime for all periods
spplot(hex_sp,c("crime"))
This map shows for each hexagon how many homicides happened. The numbers shown are for all periods.
Prepare Data Frame
#separate for each month ------------------------------------
#the previous maps was built with all periods
#now we need to separate the base for different months
crime5 = crime4 #rename to manipulate
crime5$Data = as.Date(crime5$Data,format="%d/%m/%Y")
summary(crime5$Data) ## Min. 1st Qu. Median Mean 3rd Qu. Max.
## "2008-01-01" "2011-02-09" "2014-03-16" "2014-01-01" "2016-12-06" "2019-08-22"#data max 22/08/2019 almost at the end of the period
#but I expect that this month the prediction will be bad
crime5$Data %>% class## [1] "Date"crime5$Data[sample(1:10000,5)]## [1] "2017-08-21" "2011-07-06" "2010-03-14" "2018-03-03" "2011-12-27"#create only month data
crime5$Date2 = paste0(substr(as.character(crime5$Data),1,7),"-01")
crime5$Date2[sample(1:10000,5)] %>% kable| x |
|---|
| 2013-08-01 |
| 2015-10-01 |
| 2018-09-01 |
| 2011-04-01 |
| 2014-01-01 |
crime5$Date2 %>% class## [1] "character"crime5$Date2 = as.Date(crime5$Date2,"%Y-%m-%d")
crime5$Date2 %>% class## [1] "Date"kable(crime5$Date2[sample(1:10000,5)])| x |
|---|
| 2018-01-01 |
| 2018-01-01 |
| 2011-05-01 |
| 2015-02-01 |
| 2011-01-01 |
#there are specially few homicides in august 2019, because it is not the full month
crime_c1 = as.numeric(table(crime5$Date2))
timeid = unique(crime5$Date2) #take only unique dates
timeid[1:5] ## [1] "2008-01-01" "2008-02-01" "2008-03-01" "2008-04-01" "2008-05-01"length(timeid) #there are 140 months## [1] 140crime_ts = xts(crime_c1, order.by=timeid)jpeg("time_crime.jpg", width = 800, height = 300)
plot(crime_ts)
dev.off() ## png
## 2#idea: build a spatial VAR model with time series of crimes
#need to build a data frame with each line a hexagon-region for each month
#create neighboors link
nb1 = poly2nb(hex_sp, queen=T)
#connections between hexagon-regions
plot(nb1,coordinates(hex_sp))
plot(hex_sp, add=TRUE)
This map shows the hexagon polygons and the links. Those links are used to determine the neighborhood and spatially lagged values.
jpeg("links.jpg", width = 1000, height = 1000)
plot(nb1,coordinates(hex_sp))
plot(hex_sp, add=TRUE)
dev.off() ## png
## 2W = nb2listw(nb1, glist=NULL, style ="W") #normlize weights to 1#create list with spatial data frames, with each month
df_crime = data.frame() #general data frame for models
list_poly = list() #list to input spatial polygon data frame
len_hex = length(hex_grid) #variable that contain length of spatial dataframeii = 1 #just a test
for (ii in 1:length(timeid)) {
var1 = crime5$Date2 == timeid[ii] #obs of interest
crime6 = crime5[var1,]
res = over(hex_grid, crime6[,24] , fn=sum)
res$count[is.na(res$count)] = 0
hex_sp = SpatialPolygonsDataFrame(hex_grid, data=data.frame(ID=1:len_hex), match.ID = FALSE)
hex_sp$crime = res$count
hex_sp$period = ii
hex_sp$date = timeid[ii]
hex_sp$lag_crime = lag.listw(W,hex_sp$crime)
hex_sp$long = coordinates(hex_sp)[,1]
hex_sp$lat = coordinates(hex_sp)[,2]
list_poly[[ii]] = hex_sp #input the spatial polygons in the list
names(list_poly)[ii] = ii #input name in the list for the spatial polygon
df_crime = rbind(df_crime,hex_sp@data)
#print(ii) #just to see progress
}total1=length(timeid)*len_hex; total1## [1] 15260df_crime %>% dim## [1] 15260 7#they should match
#df_crime %>% View()
df_crime[sample(1:total1,5),] %>% kable| ID | crime | period | date | lag_crime | long | lat | |
|---|---|---|---|---|---|---|---|
| 9642 | 50 | 0 | 89 | 2015-05-01 | 0.5000000 | -53.63526 | -27.05875 |
| 13961 | 9 | 0 | 129 | 2018-09-01 | 0.0000000 | -49.88526 | -28.35779 |
| 9982 | 63 | 0 | 92 | 2015-08-01 | 0.1666667 | -49.73526 | -27.05875 |
| 12499 | 73 | 0 | 115 | 2017-07-01 | 0.5000000 | -51.68526 | -26.79894 |
| 2980 | 37 | 8 | 28 | 2010-04-01 | 1.3333333 | -48.53526 | -27.57836 |
missmap(df_crime)
This is a Moran’s I test for homicides for one month. Note the p-value, there is spatial correlation in homicide.
#moran's I
var2 = list_poly[[1]]$crime
moran.test(var2,W)##
## Moran I test under randomisation
##
## data: var2
## weights: W
##
## Moran I statistic standard deviate = 2.4378, p-value = 0.007388
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.118433931 -0.009259259 0.002743691#there is spatial correlation in homicideThe following function creates a ggplot for nº of homicides and the spatial lag of homicides.
#plot graph with vec1 and spatial lag of vec1
scatLag = function(vec1,W){
vec1lag = lag.listw(W,vec1)
h2.df = as.data.frame(cbind(vec1,vec1lag))
plot = #scatterplot x=vec, y=lagged vec1
ggplot(data=h2.df, aes(y=vec1lag, x=vec1)) +
geom_point(size=.5,alpha=0.5)+
coord_equal()+
ggtitle("scatterplot x=homicide,y=lagged homicide")+
xlab("homicide") + ylab("lagged homicide")
plot
}
scatLag(var2,W)
More maps and GIFs
Here we plot how much crime there is and plot them in a GIF. This GIF shows how much crime changes in time.
# maps ang gifs -----------------------------------------
#idea: build maps of how much crime is there each regions
#build maps for each month then plot them in a gif
hex_fortify = fortify(hex_grid)
path_images = "/home/william/Dropbox/working/projects_git/Crime-SC/tests william/images/"The following loop creates all images for GIF.
ii = 1
seq1 = seq(1,length(timeid), by=6)
for (ii in seq1) {
spdf = list_poly[[ii]][,"crime"]
spdf$piece = 1:len_hex
crime_points = data.frame(coordinates(spdf), as.data.frame(spdf$crime))
names(crime_points) = c('lon', 'lat', 'crime')
crime_hex = ggmap(map, ylab = "Latitude", xlab = "Longitude") +
geom_polygon(data=hex_fortify, color='black', alpha=0, size=.1,
aes(y=lat, x=long, group=group, fill=''))+
geom_point(data=crime_points, aes(x=lon, y=lat, color=crime), size=6, alpha=.9)+
scale_colour_distiller(palette = "Spectral", limits=c(0,20))+
ggtitle(paste('crime', timeid[ii]))
ggsave(paste0(path_images,"crime_",1000+ii,".jpg"), width = 20, height = 20, units = "cm",dpi = 70)
#print(ii)
}This code chunck creates the GIF.
getwd()## [1] "/home/william/Dropbox/working/projects_git/Crime-SC/tests william"list.files(path = path_images, pattern = "*.jpg", full.names = T) %>%
map(image_read) %>% # reads each path file
image_join() %>% # joins image
image_animate(fps=1) %>% # animates, can opt for number of loops
image_write("crime_gif.gif") # write to current dirThis GIF show how number of homicides change in each hexagon. The period are changing in 6 months.
#GIF with raster ----------------------------------------------------
#build another type of GIF now with IDW function and raster format
#create raster to plot idw
pol = as(munic_dissolve,"SpatialPolygons")
grid = spsample(pol, type = 'regular', n = 90000)
gridded(grid) = TRUE
plot(grid)
We can’t see in this map. But there are 900000 squares in it. This is used in making a IDW (inverse distance weight) interpolation to create continuous maps.
grid %>% class## [1] "SpatialPixels"
## attr(,"package")
## [1] "sp"#this procedure is important to make the raster map
#idea: make all the analysis here using a square grid instead hexagonalpath_images3 = "/home/william/Dropbox/working/projects_git/Crime-SC/tests william/images3"
ii = 1
seq1 = seq(1,length(timeid), by=6)
for (ii in seq1) {
spdf = list_poly[[ii]][,"crime"]
#build raster with idw
v = as.numeric(spdf$crime)
df1 = data.frame(v)
coordinates(df1) = coordinates(spdf)
proj4string(df1) = proj4string(spdf)
names(df1@data) = c("var1")
idw = idw(var1 ~ 1, df1, grid)
idw.output = as.data.frame(idw)
names(idw.output)[1:3] = c("long", "lat", "crime")
spg = idw.output
coordinates(spg) = ~ long + lat
gridded(spg) = TRUE
raster_idw = raster(spg)
projection(raster_idw) = projection(spdf)
# map with munic borders
tm1 = tm_shape(raster_idw) +
tm_raster("crime", style = "fixed", breaks = c(0,0.5,1,2,3,4,5,6),
n = 4, palette = "Spectral", legend.show = TRUE)+
tm_shape(munic)+
tm_borders(alpha=.2, col = "black", lwd=.1)+
tm_layout(legend.text.size = .5,
legend.position = c("left","bottom"),
title=timeid[ii],
title.size=.8,
title.position = c(0.25,.3))
tmap_save(tm = tm1,
filename = paste0(path_images3,"/crime_",1000+ii,".jpg"),
width = 3, height = 3)
#print(ii) #just to follow progress
}#to build GIF for munic raster
getwd()
list.files(path = path_images3, pattern = "*.jpg", full.names = T) %>%
map(image_read) %>% # reads each path file
image_join() %>% # joins image
image_animate(fps=1) %>% # animates, can opt for number of loops
image_write("crime_gif_raster.gif") # write to current dir